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In this work we propose a criterion to finish the simulations of the Wang-Landau sampling. 
Instead of determining a final modification factor for all simulations and every lattice sizes, we 
investigate the behavior of the temperature of the peak of the specific heat during the simulations 
and finish them when this value varies bellow a given limit. As a result, different runs stop at 
different final modification factors. We apply this technique to the two-dimensional Ising model and 
a homopolymer. We verify that for the Ising model the mean order of the final modification factors 
are roughly the same for all lattice sizes, but for the homopolymer the order of the final modification 
factors increases with increasing polymer sizes. 
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I. INTRODUCTION 

Wang-Landau sampling (WLS)[H,[3 can be considered 
a well stablished Monte Carlo method, since it has been 
applied efficiently to many systems. Nevertheless the 
method is still in development and new ideas have con- 
tributed to increase the degree of efficiency and accuracy 
of the algorithm. 

The method is based on the fact that if one performs 
a random walk in energy space with a probability pro- 
portional to the reciprocal of the density of states, a 
flat histogram is generated for the energy distribution. 
Since the density of states produces huge numbers, in- 
stead of estimating g(E), the simulation is performed for 
S(E) = In g(E). At the beginning we set S(E) = for 
all energy levels. The random walk in the energy space 
runs through all energy levels from E m i n to E max with 
a probability p(E -> E') = min(exp [S(E) - S(E')], 1), 
where E and E' are the energies of the current and 
the new possible configurations. Whenever a configu- 
ration is accepted we update H(E') — > H(E') + 1 and 
S(E') -> S(E')+Fi, where F t = In/,, f = e = 2.71828... 
and /j+i = y/Jl (fi is the so-called modification factor). 
The flatness of the histogram is checked after a number 
of Monte Carlo (MC) steps and usually the histogram is 
considered flat if H(E) > 0.8(H), for all energies, where 
(H) is an average over energies. If the flatness condition 
is fulfilcd we update the modification factor to a finer one 
and reset the histogram H{E) = 0. The original version 
of WLS prescribes that simulations should be in general 
halted when / ~ 1 + 10~ 8 . Having in hand the density 
of states, one can calculate the canonical average of any 
thermodynamic variable as 



ought to update it after each Monte Carlo sweep [6|; (b) 
WLS should be carried out only up to In / = In f final de- 
fined by the canonical averages during the simulations; 
and (c) the microcanonical averages should not be ac- 
cumulated before In / < lnf m i cro defined by the micro- 
canonical averages during the simulation. The adoption 
of these easily implementable changes leads to more ac- 
curate results and saves computational time. They inves- 
tigated the behavior of the maxima of the specific heat 



C=({E-(E)j 2 )/T 2 
and the susceptibility 

X = L 2 ((m-(m)) 2 )/T, 



(2) 



(3) 



where E is the energy of the configurations and m is the 
corresponding magnetization per spin, during the WLS 
for the Ising model on a square lattice. They observed 
(as in [lOI) that a considerable part of the conventional 
Wang-Landau simulation is not very useful because the 
error saturates. In order to define ffi na i to a given model 
one should take a representative size (L = 32 for the 2D 
Ising model and N = 50 for the homopolymer) and find 
out when the corresponding canonical averages obtained 
from a few independent runs would come to steady val- 
ues. In that study they concluded that ffinai should be 
/13 and /is for the 2D Ising model and the homopolymer, 
respectively. 

In this work we propose a criterion for finishing the 
simulations which turns the choice of ffinai automatic 
for each independent run. As a result the simulated data 
become more accurate and one has no need to find out 
f final before initiating the simulations. 



(X) T = 



E E (x)Eg(E)e-? E 



(1) 



where (X)e is the microcanonical average accumulated 
during the simulations. 

Recent works have demonstrated that (a) instead 
of updating the density of states after every spin-flip, one 



II. A CRITERION FOR FINISHING THE 
SIMULATIONS 

From time to time during the WLS the random walker 
pauses the simulation in order to check the histogram for 
flatness. The number of MCS between two checks is not 
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relevant for the final resulting density of states, but it is 
interesting to relate it to the number of energy levels Ne 
of the system. In this work we used IONe for the 2D 
Ising model and lOOiVj; for the homopolymer. 

Applying WLS to the two-dimensional Ising model, be- 
ginning from fn, we calculate the temperature of the 
peak of the specific heat defined in Eq.Q using the 
current g(E) and from this time on this mean value is 
updated whenever the histogram is checked for flatness. 
When the histogram is considered flat, we save the value 
of the temperature of the peak of the specific heat T c (0). 
We then update the modification factor fi + i = \ffi and 
reset the histogram H(E) = 0. During the simulations 
with this new modification factor we continue calculat- 
ing the temperature of the peak of the specific heat T c (t) 
whenever we check the histogram for flatness and we also 
calculate the checking parameter 



\T c (t)-T c (0)\. 



(4) 



The idea of the proposed criterion is that if e remains less 
then a determined value limit during this modification 
factor, then we save the density of states and the micro- 
canonical averages and stop the simulations. In order to 
define which would be the ideal value for limit we first 
of all tested the range of variation of the temperatures. 
Defining T m i„ and T max as the minimum and maximum 
temperatures for N runs independent runs for a simula- 
tion up to /ig, we can define the mean distance between 
two temperatures as 



A = (T„ 



Tmin) i runs l) ■ 



(5) 



For the Ising model we obtained A = 6.0 x 10 4 for 
L = 32 with N runs = 24 and A = 3.7 x 10~ 4 for L = 80 
with N runs = 12- Since most of the values of temperature 
fall close to the center of the range, limit should be at 
least of the order of 2.0 x 10 -4 . Nevertheless, as we show 
in Fig. ED 



limit = 10 



(6) 



is a more stringent, but a safer parameter. Top of Fig. Q] 
shows the evolution of the temperature of the peak of the 
specific heat beginning from /n calculated for L = 80 as 
a function of the Monte Carlo sweeps @ using the 80% 
flatness criterion. The dots show where the modification 
factor was updated. The lower panel shows the evolution 
of the checking parameter e during the same simulation. 
In this example one can see that if limit = 2x 10~ 4 was 
adopted the simulation would be finished at /14, before 
the effective stabilization of the output. Using limit = 
10~ 4 the simulation proceeds up to /i6- The dotted line 
shows what would be the results if the simulation would 
be continued up to /19. 

This criterion for finishing the WLS has two main ad- 
vantages. First of all it is not necessary to determine 
f final for a representative size, as prescripted in refer- 
ence 3 since it is defined automatically for each inde- 
pendent run. The second convenience is that different 
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FIG. 1. (color online). Upper panel: Evolution of the temper- 
ature of the extremum of the specific heat during the WLS, 
beginning from /n for a single run. The dots show where the 
modification factor was updated. Lower panel: Evolution of 
the checking parameter e during the same simulation. 



runs can proceed up to different final modification fac- 
tors, depending on the evolution of the simulation. 



III. RESULTS 

A. 2D Ising model 

In order to test the proposed criterion, we performed 
Wang-Landau simulations of the two-dimensional Ising 
model with L = 32,36,40,44,48,52,56,64,72, and 80, 
taking TV = 24, 24, 20, 20, 20, 16, 16, 16, 12, and 12 inde- 
pendent runs for each size, respectively, adopting the 80% 
and the 90% flatness criteria. 



L 


f final 


f final 




80% 


90% 


32 


15.67(21) 


16.29(11) 


36 


15.25(14) 


16.71(15) 


40 


15.55(18) 


16.45(18) 


44 


15.55(23) 


16.60(17) 


48 


15.50(15) 


16.35(17) 


52 


15.06(17) 


16.25(19) 


56 


15.31(22) 


16.25(23) 


64 


15.44(22) 


16.33(22) 


80 


15.08(29) 


16.92(26) 



TABLE I. Mean order of the final modification factor for each 
simulated size of the 2D Ising model using the 80% and the 
90% flatness criterions. 
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As we pointed out before, different independent runs 
may finish at different final modification factors. For L = 
32 and 80% of flatness ffi na l ranges from /14 to fig, and 
using 90%, from /15 to /17. For L = 80 we obtained 
/14-/17 and /15-/18, respectively. 

In Table U we show the mean order of the final mod- 
ification factors for each lattice size of the simulation. 
One can see that they are roughly size independent and 
slightly more stringent then the modification factor /13 
prescribed in Ref.[3j. 

According to finite-size scaling theory [f2l - fl4j from the 
definition of the free energy one can obtain the zero field 
scaling expressions for the magnetization and the suscep- 
tibility, respectively by 

m w L-^ v M{tL^ v ), (7) 



X « I?l v X(t.l}l v ). (8) 

We see that the locations of the maxima of these func- 
tions scale asymptotically as 

T e (L) « T c + a q L~ l ' v , (9) 

where a q is a quantity-dependent constant, allowing then 
the determination of T c . 

Using these scaling functions and assuming v = 1, we 
estimated the critical temperature and the critical expo- 
nents /3 and 7. 



Case 




P 


7 


Exact 


2.2691853... 


0.125 


1.75 


Ref.[3] 


80% 


2.26934(23) 


0.1270(16) 


1.7631(27) 


90% 


2.26916(12) 


0.12494(68) 


1.7555(32) 


Present work 


80% 


2.26894(25) 


0.1239(13) 


1.7562(33) 


90% 


2.269006(77) 


0.12464(45) 


1.7572(21) 



TABLE II. Finite size scaling results for the critical temper- 
ature and the critical exponents fi and 7. 



In Table HI! we show our results for the critical temper- 
ature and the exponents /3 and 7 and compare them with 
the exact values and those of Ref. [|| . One can see that 
the accuracies for the 90% flatness criterion are roughly 
of the same order in both works, but for the 80% flat- 
ness case the adoption of the criterion for finishing the 
simulations proposed in this paper has lead to improved 
results. This improvement is welcome because the 90% 
flatness criterion is in general difficult to apply to other 
systems [l5l - fl8j . resulting sometimes in non-convergence 
or even more inaccurate values. 
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FIG. 2. (color online). Upper panel: Evolution of the temper- 
ature of the extremum of the specific heat of the self-avoiding 
homopolymer of L = 100 during the WLS, beginning from 
/16 for a single run. The dots show where the modification 
factor was updated. Lower panel: Evolution of the checking 
parameter e during the same simulation. 



B. Homopolymer 



N 


ffinal 


N 


ffinal 




80% 




80% 


40 


18.70(30) 


110 


22.30(70) 


50 


19.60(40) 


110 


23.10(69) 


60 


20.00(49) 


120 


23.50(92) 


70 


21.50(40) 


130 


25.00(30) 


80 


21.50(45) 


140 


25.30(45) 


90 


21.22(52) 


150 


24.30(45) 



TABLE III. Mean order of the final modification factor for 
each simulated size of homopolymers using the 80% flatness 
criterion. 

As a new application of the criterion, we consider a 
homopolymer consisting of N monomers which may as- 
sume any self avoiding walk (SAW) configuration on a 
two-dimensional lattice. Assuming that the polymer is in 
a bad solvent, there is an effective monomer-monomer at- 
traction in addition to the self-avoidance constraint rep- 
resenting excluded volume. For every pair of non-bonded 
nearest-neighbor monomers the energy of the polymer is 
reduced by e. The Hamiltonian for the model can be 
written as 

n = -e^2a i a j , (10) 

<i,j> 
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where a = 1(0) if the site i is occupied(vacant), and 
the sum is over nearest- neighbor pairs fl9|. (The sum is 
understood to exclude pairs of bonded segments along 
the chain.) We used the so-called reptation or "slith- 
ering snake" move which consists of randomly adding 
a monomer from one end of the chain and removing a 
monomer from the other end, mantaining the size of the 
polymer constant. We define one Monte Carlo step as N 
attempted moves. 



We performed simulations for polymers of sizes N = 
40, 50, 60, 70, 150 taking 10 independent runs for each 
size and finishing the process when the condition 




FIG. 3. (color online). Specific heat per monomer as a func- 
tion of temperature for N — 70,90, 110, 130 and 150, which 
yield ground state energies E g = —54,-72,-90,-108 and 
— 126, respectively. 



\T c (t) 



< 10" 



(11) 



was satisfied in the course of the whole simulation of a 
modification factor, where again, T c (0) is the last tem- 
perature of the peak of the specific heat in the previous 
modification factor. 



In Fig. [3] wc show the specific heat per monomer 
of the homopolymer for several chain sizes. The sim- 
ulations were started from constructed ground state 
configurations [24| . 

In both examples of models discussed here we have ob- 
tained good results by adopting limit = 10~ 4 , but when 
applying the method to more elaborate models, such as 
the HP model of protein fol ding[20l l2l| or continuous (off- 
lattice) models of polymers 22, [23[ one should always ob- 
serve the evolution of the canonical average for each new 
case, not excluding the possibility of adopting a more 
stringent parameter, such as limit = 5 x 10~ 5 , if the 
convergence is shown to be more unstable. 



IV. CONCLUSIONS 



In Fig[2] we show the evolution of the temperature of 
the extremum of the specific heat of a polymer of N = 
100, beginning from /ig. In the lower panel the checking 
parameter e of the same simulation is displayed. The 
dotted lines are the continuation of the simulation for 
one extra modification factor. Now f final varies from 
fi7 1° /20 f° r L = 40 and from / 2 2 to /26 for L = 150, 
for example. 



In Table IIIII we show the mean order of the final mod- 
ification factor for each polymer size. One can see that 
unlike the two-dimensional Ising model, in this case the 
order of the final modification factor increases with in- 
creasing polymer sizes. Adopting the proposed criterion 
for halting the simulations garantees that each particular 
run proceeds up to the real stabilization of the results. 



We proposed a criterion to finish the simulations of 
the Wang-Landau sampling. Instead of determining a fi- 
nal modification factor f final for all simulations and ev- 
ery lattice sizes, the behavior of the temperature of the 
peak of the specific heat is checked during the simulations 
and the process is halted when this value varies bellow 
a given limit during the whole simulation of a modifi- 
cation factor. As a result, different runs stop at differ- 
ent final modification factors. We applied this technique 
to the two-dimensional Ising model and a homopolymer 
and found that for the Ising model the mean final order 
of the modification factors are roughly the same for all 
lattice sizes, but for the homopolymer the final order of 
the modification factor increases with increasing polymer 
sizes. 



V. ACKNOWLEDGMENT 

This work was supported by FUNAPE-UFG. 



5 



[1] F. Wang, D. P. Landau, Phys. Rev. Lett. 86, 2050 (2001). 
[2] F. Wang and D. P. Landau, Phys. Rev. E. 64, 056101 
(2001). 

[3] A. A. Caparica and A. G. Cunha-Netto, Phys. Rev. E 

85, 04702 (2012). 
[4] L. S. Ferreira and A. A. Caparica, Int. J. Mod. Phys. C 

23, 1240012 (2012). 
[5] L. S.Ferreira, A. A. Caparica, M. A. Neto and M. D. 

Galiceanu, J. Stat. Mech. P10028 (2012). 
[6] A Monte Carlo sweep consists of L 2 spin-flip trials in the 

2D Ising model or N monomer moves in the homopoly- 

mer. 

[7] C. Zhou and J. Su , Phys. Rev. E 78, 046705 (2008). 
[8] R. E. Belardinelli and V. D. Pereyra, Phys. Rev. E 75, 
046701 (2007). 

[9] R. E. Belardinelli, S. Manzi, and V. D. Pereyra, Phys. 

Rev. E 78, 067701 (2008). 
[10] R. E. Belardinelli and V. D. Pereyra, J. Chem. Phys. 

127, 184105 (2007). 
[11] A. D. Swetnam and M. P. Allen, J. Comput. Chem. n/a. 

doi: 10. 1002/jcc. 21660. 
[12] M. E. Fisher, in Critical Phenomena, edited by M. S. 

Green (Academic, New York, 1971). 
[13] M. E. Fisher and M. N. Barber, Phys. Rev. Lett. 28, 

1516 (1972). 



[14] in Phase Transitions and Critical Phenomena, edited by 
C. Domb and J. L. Lebowitz (Academic, New York, 
1974), Vol. 8. 

[15] C. J. DaSilva, A. A. Caparica and J. A. Plascak, Phys. 

Rev. E 73, 036702 (2006). 
[16] C. J. DaSilva, A. G. Cunha-Netto, A. A. Caparica and 

R. Dickman, Braz. J. Phys. 36, no. 3A 619 (2006). 
[17] A. G. Cunha-Netto, R. Dickman and A. A. Caparica, 

Comput. Phys. Comm. 180, 583 (2009). 
[18] D. T. Seaton, T. Wiist, and D. P. Landau, Phys. Rev. E, 

81, 011802 (2010). 
[19] R. Dickman, J. Chem. Phys. 96, 1516 (1992). 
[20] T. Wiist and D. P. Landau, Comput. Phys. Comm. 179, 

124 (2008). 

[21] T. Wiist and D. P. Landau, J. Chem. Phys. 137, 064903 
(2012). 

[22] T. Wiist and D. P. Landau, Phys. Rev. Lett. 102, 178101 
(2009). 

[23] M. P. Taylor, W. Paul and K. Binder, J. Chem. Phys. 

131, 114907 (2009). 
[24] C. J. DaSilva, A. G. Cunha-Netto, R. Dickman and A. 

A. Caparica, Comput. Phys. Comm. 180, 590 (2009). 



